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We investigate the resonantly enhanced tunneling dynamics of ultracold bosons loaded on a tilted 
1-D optical lattice, which can be used to simulate a chain of Ising spins and associated quantum 
phase transitions. The center of mass motion after a sudden tilt both at commensurate and incom- 
mensurate fillings is obtained via analytic, time-dependent exact diagonalization and density matrix 
renormalization group methods (adaptive t-DMRG). We identify a maximum in the amplitude of 
the center of mass oscillations at the quantum critical point of the effective spin system. For the 
dynamics of incommensurate systems, which cannot be mapped to a spin model, we develop an 
analytical approach in which the time evolution is obtained by projecting onto resonant families of 
small clusters. We compare the results of this approach at low fillings to the exact time evolution 
and find good agreement even at filling factors as large as 2/3. Using this projection onto small 
clusters, we propose a controllable transport scheme applicable in the context of Atomtronic devices 
on optical lattices ('slinky scheme'). 
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I. INTRODUCTION 

Ultracold atoms on optical lattices provide a well con- 
trolled system for the study of out-of-equilibrium dynam- 
ics, due to the experimental ability to tune their mi- 
croscopic parameters, and to even do this dynamically 
during the course of a single experiment [1]. In recent 
years it has been proposed to use these clean and tun- 
able systems for the realization of ultracold atom-based 
analogs of semiconductor electronic devices (Atomtron- 
ics) [2, 3]. Recent developments in single-site addressing 
[4-6] and in the ability to engineer optical lattices of ar- 
bitrary geometry are important steps towards achieving 
this goal. At the heart of Atomtronics is quantum trans- 
port. In condensed-matter-based electronic devices the 
simplest way to obtain transport is to apply an electric 
field. However, for pure systems in a periodic potential, a 
constant force generally leads to Bloch oscillations (BO) 
[7]. Except from special cases such as clean semiconduc- 
tor superlattices [8], BO are not an issue in the complex 
solid state environment since they are quickly damped 
due to strong dissipative processes, e.g., defects, impu- 
rity scattering, and band tunneling. Ultracold gases are 
on the contrary almost perfectly decoupled from the en- 
vironment and therefore BO dominate the dynamics for 
long times [9, 10]. These considerations pose the ques- 
tion as to how one may induce transport in an optical 
lattice while keeping full control over the atomic quan- 
tum states, as required in Atomtronics. 

Previous studies aiming to enhance transport on opti- 
cal lattices have mainly focused on the regime of weak in- 
teractions in which the particles are delocalized along the 
system [11-14]. One reason for this is that strong interac- 
tions tend to localize particles and thus inhibit transport. 
Nevertheless, quantum transport in the strongly interact- 
ing regime has attracted some interest [15-18] with a fo- 



cus on gaining further insights into transport of quantum 
information through a lattice system and quantum com- 
munication in network systems. Here, however, we ad- 
dress directly the problem of enhancing transport of par- 
ticles through a lattice system by considering the strongly 
interacting regime of the repulsive Bose-Hubbard model 
(BHM) with an external linear potential. We show that 
by resonantly tuning the strength of the linear poten- 
tial, i.e. by adjusting the bias between adjacent wells 
to allow tunneling within the lowest band, it is possible 
to highly control and understand, even analytically, the 
complicated many-body dynamics. We assume a deep 
enough lattice with suppressed interband tunneling due 
to a large band gap. We emphasize that what we denote 
as resonance should not be confused with the resonances 
originated by tunneling to excited Wannier-Stark lad- 
ders discussed in Ref. 19. Based on this understanding, 
we propose an approach for enhancing transport, which 
we call the 'slinky scheme' due to the peculiar nature 
of the resulting motion. Somewhat counterintuitively, it 
is the presence of strong interactions which leads to an 
enhanced particle motion. 

In addition to studying BO and transport properties, 
strongly interacting bosons on a tilted lattice at commen- 
surate fillings have been predicted to be a very fruitful 
system for the investigation of quantum magnetism in 
Ising models and quantum phase transitions [20]. The 
underlying idea is to map the doublon-hole excitations of 
a Mott insulator onto an effective spin degree of freedom 
when the applied tilt is tuned resonantly to the doublon 
interaction energy. Just recently the observation of the 
associated quantum Ising transition has been reported 
in the laboratory [21], leading to proposals for the re- 
alization of other systems, e.g., quantum dimer models 
[20, 22]. These experiments also motivate the study of 
the time evolution in the resulting spin systems. In this 
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paper, we do so by characterizing and comparing the dy- 
namics of the center of mass (CM) oscillations of both the 
BHM and the effective Ising system. This is achieved by 
numerically computing the time evolution of both models 
using time-dependent exact diagonalization (ED) [23-26] 
and adaptive time dependent density matrix renormaliza- 
tion group (adaptive t-DMRG) [27-31] techniques. The 
key result is that there is a maximum in the amplitude 
of the CM oscillation at the quantum critical point. We 
extend our studies also to incommensurate fillings where 
the mapping to an Ising spin model does not hold, and 
we develop an analytical framework to deal with the dy- 
namics in the low-density regime. The method relies on 
projecting the time evolution onto subspaces spanned by 
resonant families of states on small clusters. We com- 
pare the approximate results to t-DMRG results and find 
that the analytical approach provides a good approxima- 
tion for filling factors as large as 2/3. Using this simple 
method one can show that the CM oscillations are re- 
stricted to have an amplitude of less than one lattice 
spacing. Even though the amplitude increases with the 
filling factor and is maximal at unit filling, transport re- 
mains inhibited. To overcome this limitation, we propose 
an alternative scheme. We first pattern load atoms using 
a superlattice [32] and then stroboscopically modulate 
the amplitude of the lattice. The resulting motion of the 
atoms is reminiscent of the motion of a 'slinky' toy down 
a stairway. We expect this approach to be realizable with 
available experimental methods. 

The paper is organized as follows. In Sec. II, we intro- 
duce the model, the non-equilibrium set-up, the mech- 
anisms underlying the BO, and the methods used. In 
Sees. Ill and IV we discuss the center of mass motion at 
commensurate and incommensurate fillings, respectively. 
In Sec. V we discuss our numerical results obtained via 
time-dependent ED and adaptive t-DMRG. We discuss 
a possible signature of the critical point of the spin sys- 
tem in the amplitude of the CM oscillations in Sec. VB. 
In Sec. VI we present the transport scheme mentioned 
above. In Sec. VII, we conclude with some prospects for 
future work. 



II. MODEL, NON-EQUILIBRIUM SETUP AND 
METHODS 



A. Model and Center-of-Mass observable 

We treat the one-dimensional single-band Bose Hub- 
bard model on a tilted lattice, 



{i-j) » i 

where () denotes the sum over neighboring lattice sites, 
) is the annihilation (creation) operator for a boson at 
site i, and n,- = a\a.i, and with f2 the strength of the tilt- 
ing potential. This single-band model represents a good 



approximation in the presence of a deep lattice which pre- 
vents inter-band transitions and decay of the atoms out 
of the lattice [19]. There are several ways of realizing the 
tilted lattice in experiments. For example, it is possible 
to exploit the gravitational potential by creating a verti- 
cal optical lattice, or to detune the counter-propagating 
laser-beams forming the optical lattice. In the latter case, 
a time dependent detuning 5^{t) can lead to an acceler- 
ation in the lattice depending on the induced velocity 
v{t) = \5^{t)/2 [33], with A the wavelength of the laser. 
The most important observable for the treatment of the 
Bloch oscillations is the time evolution of the center of 
mass (CM) position 
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with N the total number of particles on the lattice with L 
sites, and Xcm is measured in units of the lattice spacing, 
i.e., the lattice spacing d in all subsequent calculations 
is set to unity. Measuring Xcm in the experiments [34- 
38] provides not only information on the CM motion and 
the associated current (icm(^) = —^{[xtH]), where h is 
equal to one throughout the paper) significant to trans- 
port properties, but can also provide insights into rele- 
vant energy scales in the system. For instance, the CM 
motion in Bloch oscillations has proven to be an accurate 
tool for metrology [10, 39, 40]. 



B. Bloch Oscillations 

In this section, we review the main properties of Bloch 
oscillations (BO) of a single particle on a tight-binding 
chain subjected to a linear tilt. In this case, the en- 
ergy levels are discrete e(n) = fin and the eigenstates are 
called Wannier-Stark states [19, 41]. In the Wannier basis 
they are given by \4>n) = Jj-n{'2.J/fl)\j), where Jm{x) 
is the Bessel function of the first kind. The nth Wannier- 
Stark state is localized in the region \ j — n\< 2 J/17. Due 
to the harmonic-oscillator like spectrum, an initial state 
centered at a;cm(0) [42] exhibits periodic CM oscillations 
with frequency Q, 
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The same behavior is also captured by a semiclassical 
picture in which the linear potential is treated as a con- 
stant force dragging the particle through the Brillouin 
zone. This treatment leads to a time-dependent value 
of k and hence to a group velocity and to corresponding 
CM oscillations of the form 

k^n^ k{t) = fco + ^t, 
, , dElklt)) , ^ 2J ^ (4) 

Vg{t) = ^^j^ ^ XcM{t) - — COS fit, 

with E{k) = — 2Jcos(fc) the dispersion of the system 
without the external potential. 
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A system of many non-interacting bosons on the lat- 
tice will show the same dynamics. Inter-particle interac- 
tions, in contrast, tend to dampen the CM motion. In 
the weakly interacting regime, interactions lead to un- 
derdamped dynamics [43, 44], while in strong fields they 
can yield a multitude of interesting phenomena, e.g.^ BO 
at interaction-induced frequencies or so-called quantum 
carpets [45, 46]. In the strongly interacting regime, the 
CM oscillations are almost completely suppressed. How- 
ever, time-dependent interactions can in principle stabi- 
lize them [47] . 

In the limit [/ — )■ oo, a hard-core constraint is realized 
which mimics the Pauli exclusion principle since double 
occupancy is suppressed. In this fermionized regime it 
is then possible to introduce a Jordan- Wigner transfor- 
— i=i '^i'^i , which reduces the Hamilto- 

nian to a tight-binding model for non-interacting spin- 
less fermions. Any local observable such as the CM 
position is identical for spinless fermions and hard-core 
bosons. If we denote by |?A*(0)) the initially populated 
single particle states with q = 1, - ■ ■ N, the evolution of 
each of them in the basis of Wannier-Stark states reads 

\r{t)) = E„ /«e-»"* |0„), with = (M^m. Us- 

ing the recurrence relation of Bessel functions (see Ap- 
pendix A) for the CM oscillations of TV hard-core bosons 
on an infinite lattice, one obtains 
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average position x = ^EgE„'^l/«l 
find that at larger fillings the motion is suppressed, as 
demonstrated in Fig. 1. The simple physical picture for 
this effect is that as the filling increases hard-core bosons 
have on average less space for free motion before they 
encounter each other. As the system approaches unit 
filling, particle transport gets fully suppressed along the 
lattice. 



C. Methods 

We apply different approaches for the treatment of 
the dynamics of the system Eq. (1). For systems small 
enough, we compute the time evolution by fully diag- 
onalizing the Hamiltonian matrix. For larger systems, 
we apply a Krylov-space approach to the time evolution 
in the framework of exact diagonalization (Krylov-ED). 
For the largest system sizes treated, we apply the Krylov- 
space variant of the adaptive time dependent density ma- 
trix renormalization group method (adaptive t-DMRG). 
In all cases, we introduce a cutoff in the local bases on 
the lattice sites and keep up to 3 bosons per site. In the 
Krylov-ED approach, we approximate the time evolution 
operator in a basis of = 10 Lanczos vectors and use 
a time step of At = 0.0005 for the BHM. The resulting 
error on the time scales treated is typically of the order 
of machine precision, and we can treat systems up to 



15 sites with 8 particles. With the adaptive t-DMRG, 
we aim for a discarded weight < 10~^ and keep up to 
m = 500 density-matrix eigenstates during the time evo- 
lution. For the BHM, we use a time step At = 0.0005 
and find at the end of the time evolutions displayed in 
the plots a discarded weight < 10^^. The effective spin 
models are much easier to treat, and we apply a time step 
of At = 0.005, resulting in a discarded weight of < 10~^ 
at the end of the time evolution for our largest system 
size of L = 50 sites. 



III. CENTER OF MASS OSCILLATIONS: 
COMMENSURATE FILLING 

In this section we study the CM oscillations at com- 
mensurate filling {N = L), introduce the concept of res- 
onant dynamics and compare it to the BO. We start the 
discussion with the pedagogical example of a double- well 
system with strong interactions U ^ J. 



A. The simplest system: a double well potential 

We discuss the basic concept and properties of a dipole 
[20] in a tilted double well. We assume that the system 
is initially in a unit-filled state |0) = |11) with one par- 
ticle in each well. In the absence of a bias, this "Mott 
insulating" state is to a good approximation the ground 
state for U ^ J. Excited states are obtained by moving 
one particle to the other well, i.e., we obtain particle-hole 
excitations with two particles on one site and no parti- 
cles in the other. When applying a tilt, the lowest lying 
excitation is the one in which the two bosons are in the 
lower potential well. This situation can be modeled by 
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FIG. 1. (Color online) Center-of-mass motion due to Bloch 
oscillations for the BHM Eq. (1) in the hard-core limit for dif- 
ferent fillings n — N/L, where L = 40 sites and Q = 40 J. The 
time evolution is obtained by diagonalizing the corresponding 
single-particle Hamiltonian and we show results for n — 1/40 
(Blue, short dashes), n = 1/4 (Red, long dashes) and n = 1/2 
(Green, dashed-dotted line). 
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FIG. 2. (Color online) CM motion of a double well system at 
resonance. The inset shows the energy levels of the different 
states as a function of detuning (U — Q.)/J. 



applied to the unit-filled ground state of the untilted po- 
tential. For a generic bias < J7, the particles remain 
localized on each site, since the hopping can overcome 
the energy cost of neither the potential nor of the on-site 
interaction (see the inset of Fig. 2). However, near reso- 
nance ([7 = 0), the unit-filled state and the dipole state 
are nearly degenerate, and so tunneling from one well to 
the other is possible. 

The time evolution of the double-well system is given 
by \il){t)) = co(i)|ll) + ci(i)|20) + C2(t)|02), and can be 
obtained analytically by an approach similar to the adi- 
abatic elimination of a non-resonantly coupled excited 
state in a three-level lambda system [48, 49]. At res- 
onance, perturbation theory shows that the population 
of the off-resonant state |02) is ^ J^/il^ and hence for 
n ^ J remains vanishingly small during the time evolu- 
tion (see Appendix B). To lowest order, the CM motion 
is described approximately by 
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(6) 



The high frequency oscillation ~ is due to the popu- 
lation of the non-resonant state separated by a large en- 
ergy difference from the other states (Fig. 2). Since the 
population of this state is strongly suppressed, setting 
C2{t) = is a good approximation. Therefore, all terms 
~ jVt} are neglected, including the high frequency os- 
cillations. 

This simple double well case illustrates that the time 
evolution of the system can be obtained in a good approx- 
imation by neglecting contributions from non-resonant 
states. This will be used throughout the rest of the pa- 
per. 



B. Generic case 

We now study the behavior of the CM oscillations as 
a function of U while keeping J and S> J fixed. Fig- 
ure 3 summarizes the main results. At ?7 = the CM 
shows BO with amplitude 4J/ri and frequency fi. For 
our choice of J and ^ J even at [/ = the BO's am- 
plitude is less than one lattice spacing. Upon increasing 
J7, first the amplitude of the BO decreases by a factor of 
more than 2 for JJ/fi = 0.5, in accord with the general 
expectation that BO are suppressed in the presence of in- 
teractions. However, further increasing U and tuning the 
system to resonance, J7 = f2, leads to enhanced CM oscil- 
lations. The period of these oscillations scales as J rather 
than fi, as in the case of BO, indicating that the mecha- 
nism underlying the dynamics is very different from the 
one of the BO [50]. Note that even though the ampli- 
tude of the CM oscillations is maximal at resonance, it 
is < 1/2 lattice spacings, regardless of the system size. 
This is due to the fact that the bosons are restricted to 
hop only between nearest neighboring sites. Since we an- 
alyze the resonance situation deep in the Mott-insulating 
regime, we require strong fields f2 ^ J. We leave out the 
resonance situation for weak i7 for future investigations. 



C. Resonant case: Effective spin model 

At resonance the dynamics is obtained by restricting to 
the set of dipole excitations that are degenerate in energy 
to zeroth order in the hopping. We display an example for 
such a resonant family of states in Fig. 4. Interestingly, 
these states have a representation in terms of an effective 
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FIG. 3. (Color online) CM motion for a BHM at commensu- 
rate filling with L = 15 lattice sites for J = 1, S7 = 40 and 
three different values of U . The results at f/ = are obtained 
by diagonalizing the single particle problem, the results at fi- 
nite \J are obtained by adaptive t-DMRG, with a discarded 
weight < 10~® at the end of the time evolution. The results 
for the non-interacting system show the frequency f2 and the 
amplitude 2,J/Q. of the BO, while at resonance the frequency 
is ~ J and the amplitude gets strongly enhanced. 
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spin language which we discuss below. In this picture, the 
resonant states map to all spin configurations, excluding 
the ones containing two adjacent | f) states. 

The equilibrium properties of this resonant family of 
states have been derived in Ref. 20. For the sake of clar- 
ity, we summarize them in this section. 

We work in the effective basis obtained by all possible 
coverings of the system by dipoles, spanned by all com- 
binations \M;k) = YifLi '^l: \^) ^ where |0) is the initial 
Mott-insulating state with one particle per site and M 
the number of dipoles in the system; the vector k contains 
the positions of the M dipoles on the system. Note that 
two hard-core constraints appear: first, by construction, 
two dipoles cannot occupy the same site, leading to an 
on-site hard-core constraint. Second, it is not possible to 
realize two dipoles on adjacent sites, since the creation 
of the second dipole would destroy the first one, leading 
to a nearest-neighbor hard-core constraint. Taking both 
constraints into account, the effective Hamiltonian [20] 

HeSMp = -JV2j2idl + di) + {U-n)Y, d-U^ (7) 

i i 

is obtained. This model can be further mapped to a 
5 = 1/2 Ising model in a longitudinal and transverse 
external magnetic field, 

iJeff = - JV25] af + {U-^) Y,{at + l)/2 

(8) 

+ A^K + l)K+i + l), 

i 

by introducing pseudospin raising and lowering operators 
via the mapping d| — (Tj"^ (erf and erf are the correspond- 
ing Pauli matrices for the spin on site i). This mapping 




FIG. 4. Connected resonant family of states for a system 
with L = A lattice sites at commensurate filling. The Mott- 
insulating state in (a) is degenerate in energy to zeroth order 
in J to the state configurations with one dipole excitation in 
(b) and the state configuration with two dipole excitations 
in (c). In the spin representation, the unit filled state in 
(a) is represented by j iii). The configurations in (b) from 
top to bottom are represented by | \ and | ||t) 

respectively. The configuration in (c) is represented by | tit) • 
Note that in the resonant family of states configurations with 
two adjacent | t) states are excluded. 



captures the on-site hard-core constraint. The nearest- 
neighbor hard-core constraint is captured by choosing 
A sufficiently large, so that two neighboring spins can- 
not simultaneously point up. The quantum critical be- 
havior of this system is governed by a single parame- 
ter, A = and possesses an Ising critical point at 
Ac ~ —1.85 at which the system undergoes a phase tran- 
sition from a paramagnet to an antiferromagnet, as dis- 
cussed in Ref. 20. 

In this mapping, the spins are located on the bonds be- 
tween two sites. In addition, the dimension of the Hilbert 
space needed to obtain the time evolution via this ef- 
fective model is reduced to dim('H) = Yl^M=a — 
F{L + 1), where N^^^ = {^~m ') is the number of states 
in a system of L sites containing M dipoles and F{L) 
denotes the L*'' Fibonacci number. This has to be con- 
trasted with the dimension of the Hilbert space of the 
original model, dim(?^) — (^^^^^)- Note that due to 
the hard-core constraints the maximum number of pos- 
sible dipoles in a system is the floor function applied to 
half the system size, \L/2\. 



1. Dynamics in a basis of symmetric states 

One of the goals of the present paper is to identify 
possibly simple analytical approaches to the dynamics 
of the tilted Mott insulator. In the preceding section, 
we have already achieved a substantial simplification by 
mapping the complicated bosonic system to a relatively 
simple effective spin model. However, due to the hard- 
core constraint, the model in Eq. (8) is a many-body 
model with competing interactions, so that obtaining the 
dynamics is still a challenging task. In this section, we 
will discuss how the short time dynamics can be obtained 
analytically by introducing further approximations, and 
how the validity of this approximation can be enlarged 
by properly accounting for the hard-core constraint. 

The simplest approach is to completely neglect the 
A-term in Eq. (8). This term is essentially excluding 
the basis states with two adjacent | f)- For an initial 
state in which all spins point downwards, for times short 
enough, these basis states will remain unpopulated and 
neglecting the hard-core constraint should be a good ap- 
proximation. This treatment should also be a better 
approximation for larger detunings (|J7 — r2| 3> J). In 
this case, dipoles become more energetically costly to 
populate and the system remains in the no dipole or a 
single dipole manifold in which the hard-core constraint 
is not relevant. This is a favorable situation: without 
the A-term, we are dealing with a non-interacting sys- 
tem which can be treated exactly. The dynamics in this 
case is obtained in a basis of Dicke states [51], which 
are the set of spin states with maximum total spin. In 
this approximation, the Hamiltonian can be rewritten as 
« -2 Jy2 5- + iU-n) 5f„,,i with collective spin 
operators S^^^^i = | J^i '^i' ■ The dynamics is then de- 
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scribed as a rotation of the Bloch vector of the collective 
spin state manifold, S'totai = {L—l)/2. It can be obtained 
by solving the equations of motion for the total spin com- 
ponents with initial condition (^^^^^^[(O)) = — (i — l)/2. 
Defining wq = \/'B>,P + {U — fi)^, one then obtains 
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Using the mapping discussed later in Eq. (19), the CM 
motion is, 
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When comparing the result of Eq. (10) to the exact re- 
sults of Fig. 6, we find good agreement up to times 
Jt K, 0.3. As discussed above, the larger the detuning 
from resonance, the better the qualitative agreement. 

This approximation hence leads to a closed expression 
for the CM motion. However, neglecting the A-term is 
a very rough approximation. A more accurate treatment 
can be obtained by excluding all states with two adja- 
cent It)- In this way, we account for the hard-core con- 
straint but stay within the collective spin manifold. This 
basis we refer to as the set of symmetric states. For 
a two site system (one spin), the Dicke and symmetric 
states are trivially the same and the time-evolved state 
is spanned by either set, since the hard-core constraint 
does not come into play. For a three site system (two 
spins), the dynamics obtained in the symmetric states is 
exact. However, in the general case, the hard-core con- 
straint leads to a time-dependent phase for the different 
spin configurations, and the system decays out of the 
subspace of symmetric states. Despite of this, for times 
short enough, most of the weight of the many-body wave- 
function is on the symmetric manifold and one can treat 
the dynamics on this time scale to a good approximation. 
This is done numerically. As expected, we find that the 
result is a better approximation than the treament in the 
Dicke states, the dynamics being comparable to the ex- 
act treament for longer times. Nevertheless, the decay 
out of the symmetric manifold grows, in lowest order, 
as Pasym ^ n!^^^ J^t^ and therefore the time interval in 
which this approximation is valid shrinks with increasing 
system size. 



IV. CENTER OF MASS OSCILLATIONS: 
INCOMMENSURATE FILLING 

We now develop an approximate analytical solution to 
the dynamics of the system at incommensurate fillings, 
A'' < L, at which the spin model is not valid. We begin 



with an illustrative example by considering a lattice of 
six sites initially prepared in the Fock state |110111). 
Suppose we time evolve this state at resonance, C/ = il, 
with [/, f7 ^ J. As we discussed before, the particles 
explore configurations degenerate in energy to the initial 
state in zcroth order in J. The family of states which 
needs to be considered during the time evolution is given 
by the states 
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The time evolution of the system in the initial state \ai) 
is then 
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This results in the CM motion 
1 r 



•^cm(^) 



10 



34 + cos (2\/2Jt) -f cos(4Jt) . (13) 



1.5 • 



ci 1.0 



• 111) 
■ nil) 

♦mil) 

Allllll) 



0.5 



O.OL 



: t t * I : I I 



10 12 14 16 18 20 
Number of Sites L 

FIG. 5. (Color online) Weighted sums of the coefficients of the 
ground state of the BHM, Eq. (1) in the hard-core limit U — > 
oo for fixed A^ = 8 as a function of system size L. We consider 
configurations that possess n = 2, 3, 4, and 5 adjacently 
occupied sites (clusters), and weight the coefficients by the 
number of occurrences D^"l of the cluster which are obtained 

{k) 

numerically as described in Appendix C. 
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The form of Xcui{t) indicates that the dynamics is gov- 
erned by the two subspaces spanned by the resonant fam- 
ihes of the |11) and |111) states. Projected onto those 
subspaces the effective Hamiltonian becomes 



H' ^ pDHeffPa -V2J (di + di + + h.c. 
= Hi + H4R, 



(14) 



where the operator Pjj projects onto the resonant fami- 
hes. Due to the fact that [Hi,H4^5] = 0, the time evo- 
lution takes place in independently evolving subspaces, 



= e-'^i*|ll) (g) |0) (ge-*^*'=|lll), 



(15) 



and the total CM dynamics reduces to the direct sum of 
the CM evolution in each of the clusters: 



1 






(1 


N 




1 






1 






1 




10 





5 + cos \2V2jVj + 29 + cos (4 Jt) 



(16) 



We have thus demonstrated that the resonant dynamics 
of a larger system can be treated by considering the time 
evolution of smaller decoupled systems. In order to use 
this approach for a system of arbitrary size, we have to 
consider configurations of particles which are sparse and 
spread out over the lattice. For simplicity we will assume 
the atoms are prepared in the ground state in the absence 
of a tilt , i.e ri = 0, we assume open boundary conditions 
and that at time i = the system is suddenly tilted close 
to resonance, U ^ il ^ J. 

To be more quantitative, we rewrite the initial wave 
function in a more tractable, albeit approximate way: 



{k} 



n 



(17) 



{kr^}:{k,} 



where the summation extends over all permutations of 
the positions {ki},{kj}, {km}, and {fc;}; the coefficients 
are derived in Appendix C. |0) is the vacuum state 



and the operators aj.^ , dj,. , ^1, , and ■^l^^ create config- 



urations |1), |11), |111), and |1111), respectively, with a 
given weight determined by the coefhcient c^j:y We use 

the convention that the components of k denote the par- 
ticle positions and the subscripts attached to the creation 
operators denote the location of the leftmost site of the 



cluster. We account for subspaces up to only |1111) be- 
cause the configurations containing larger clusters can be 
presumed to carry negligible weight (sparse filling condi- 
tion, see Fig. 5) . This can be understood in the hard-core 
regime where the atoms in the initial many-body ground 
state want to spread out symmetrically with respect to 
the center of the lattice in order to maximize their kinetic 
energy, thus avoiding the energetically forbidden double 
occupancy. Hence, at low filling factors, the configura- 
tions that most significantly contribute to the ground 
state are those clusters with the lowest number of con- 
tiguous occupied sites (see Fig. 5). During the course 
of the time evolution, each cluster evolves independently 
within its resonant manifold, allowing us to treat the time 
evolution of the full system by computing the time evo- 
lution of each of the small clusters and their associated 
resonant families. In addition, the configurations |1) do 
not contribute to the dynamics since the resonance con- 
dition is significant only if at least two adjacent sites are 
occupied. 

We denote the number of n-particle clusters of a basis 
state which has particles at positions k by dI"^ . Then 
we obtain for the CM motion in this approximation 



2Nx{t) ^ |c 
{k} 



- |2/nl"> 
{k}\ \^{k} 



cos 2V2Jt 



-D\^}^^ cos (4 Jt) 

r)|iiii> 
{fc} 
34 



1 - V 17 cos 2\/ 5 - 4V17 Jt 



1 + VU] cos 2\/ 5 + AVUJt 



64 cos V 5 + 4VT7Jt cos V 5 - 4^1 7 Jt 



(18) 

where we denote by {k} all permutations of the posi- 
tions of N particles. This expression is one of the main 
results of the paper. At low enough fillings, it represents 
a good approximation to the time evolution of systems 
of arbitrary size and number of particles and predicts 
the dynamics of tilted Mott insulators at resonance away 
from the commensurable case. In the next section, we 
will test the accuracy of this approach. 



V. COMPARISON OF NUMERICAL RESULTS 
TO THE EFFECTIVE SPIN MODEL AND THE 
ANALYTICAL TREATMENT 

In this section, we compare the time evolution obtained 
via the effective spin model Eq. (8) to the exact one gov- 
erned by the BHM Eq. (1). In addition, we analyze the 
dynamics of the effective spin model Eq. (8) as a func- 
tion of A and find that the amplitudes of the BO possess 
a maximum at A,.. 
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FIG. 6. (Color online) Numerical results (Krylov-ED) for the 
CM motion in a system with L = 12 sites for three different 
values of A = ([/ — f2)/J, where J = 1 and Q, = 40. The solid 
lines show the evolution of the effective spin model, while the 
dashed lines show the exact evolution of the BHM. As the 
parameters are tuned off resonance, the amplitude decreases, 
indicating less participation of higher dipole excitations. 



A. Accuracy of the effective spin dynamics at 
commensurate filling 

In order to compare the dynamics of the effective spin 
model to the one of the BHM at commensurate filling, 
we use the mapping 



left boundary site: ni {af + 3)/2 
right boundary site: n^^/ — > (1 — (7£)/2 (19) 
bulk: n, ~> {al - ct^_i + 2) /2 

where L is the size of the spin system and L' = L + 1 the 
size of the bosonic system, which has one lattice site more 
than the spin system, due to the fact that the dipoles are 
located on the bonds. This mapping is obtained by com- 
paring individual site occupations in the presence or ab- 
sence of a dipole; the different mapping at the boundary 
sites and in the bulk is due to the open boundary condi- 
tions (OBC). The difference between the mapping at the 
left and the right boundaries is due to the fact that the 
site at the highest potential (right boundary) will pos- 
sess either zero or one particle, but the site at the lowest 
potential (left boundary) has either one or two particles. 

The dynamics of the BHM is obtained by evolving an 
initial Mott insulating ground state with one particle per 
site. In the spin picture, this is equivalent to an initial 
state with all spins pointing downwards. In Fig. 6 we 
show the time evolution of the BHM and the one of the 
effective spin model at A = 0, A = —1, and A — —5. We 
find that the effective spin model reproduces the dynam- 
ics of the BHM even when detuning not too far off reso- 
nance. The evolution of the systems is essentially iden- 
tical on short time scales (Jt < 2.5 for A = 0, < 0.5 
for A = —1, and Jt < 0.3 for A — —5) but then differs 
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FIG. 7. (Color online) Amplitudes of the CM motion as ob- 
tained in the effective spin model and as defined in the inset 
in (a), which shows the evolution at resonance. The results 
shown are for systems with L = 20, 30 and 50 lattice sites 
and are obtained using the adaptive t-DMRG. The ampli- 
tudes are obtained from interpolating the discrete data points 
of the CM motion, (a) Difference Al (measured in units of 
lattice spacings) between the initial CM position and the first 
minimum as a function of A. (b) The second amplitude A2 
(measured in units of lattice spacings) as a function of A. The 
inset shows the finite size extrapolation at A = —1.85, and the 
black circles show A2{X) after finite size extrapolation. We 
estimate the error of the extrapolation to be of the order of 
the symbol size. 



increasingly in the course of the time evolution. In addi- 
tion, for small detunings from resonance, the frequency 
of the oscillation is in excellent agreement to the exact 
solution, while it differs in the case of stronger detun- 
ings. Note that detuning causes the system to be rigid in 
the sense that states that contain more dipoles are more 
costly in energy. Hence, on sufficiently short times, the 
number of dipole states contributing to the dynamics is 
reduced so that the time evolution can be obtained by 
considering a smaller number of states (see Ref. 52 for a 
detailed discussion). 



B. Maximum of the CM amplitude at Ac 

Now we turn to the properties of the CM oscillations 
in the time evolution of the effective spin system Eq. (8) 
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FIG. 8. (Color online) (a) CM motion for a system of L = 15 
lattice sites at resonance at low fillings away from commen- 
surability. The solid curves show the result of the analytic 
calculation, Eq. (18). Note that the analytic expression ne- 
glects the fast BO which are caused by population of non 
resonant states, (b) Spectral analysis of the CM motion. At 
low fillings, 2 particle cluster states dominate the dynamics. 
The high frequency BO are visible as peaks at il/J = 40 
which decrease upon increasing the filling. The additional 
low frequency peaks are due to larger cluster states which be- 
come relevant at higher fillings. The results for the BHM at 
n = 2/5 and n = 8/15 are obtained via Krylov-ED, the ones 
at n = 2/3 are adaptive t-DMRG results. 



upon changing A. As before, we consider the time evo- 
lution of an initial state in which all spins are pointing 
downwards, which is equivalent to a Mott-insulator with 
one particle per site in the bosonic language. As shown 
in Ref. 20, the effective spin model possesses a quantum 
critical point at Ac ~ —1.85. We are interested in pos- 
sible signatures of this critical point in the CM motion. 
In Fig. 7 we present our adaptive t-DMRG results for 
the first oscillations Al, A2 and A3 [defined in the in- 
set of Fig. 7(a)] as a function of detuning A for systems 
with L = 20, 30 and 50 sites. Interestingly, we find a 
local maximum of A2 at A ~ —1.85, indicating that this 
quantity might indeed reveal the existence of a critical 
point in the time evolution of this system. This fea- 
ture seems to persist upon changing the system size, and 
at Ac we obtain after finite size extrapolation the value 
A2{Xc,L — )■ oo) w 0.255. This can be contrasted to the 



value of A2(-0.05, L ^- oo) ss 0.227, so that the am- 
plitude when increasing the detuning from resonance at 
A = to Ac changes by « 10%. Similarly, we find a max- 
imum of AS at A « —1.85. However, such a maximum in 
the vicinity of Ac does not show up in Al(X), a quantity 
which instead reaches its peak value at resonance A = 0. 

This puts forth the interesting possibility that one may 
be able to use the CM oscillations to identify quantum 
critical points of a generic phase transition. 



C. Accuracy of the effective cluster dynamics at 
incommensurate fillings 

In Fig. 8 we compare the CM motion as obtained 
by the original BHM for a system of L = 15 sites at 
various fillings n < 2/3 to the one obtained in the 
approximate treatment using families of small cluster 
states. At fillings n = 2/5 and n = 8/15 the plot shows 
data obtained via the Krylov-ED approach, so that the 
results are essentially exact. At n = 2/3, we present 
data obtained via adaptive t-DMRG where the discarded 
weight is > 10^^ for times Jt > 2, so that the results 
at later times might be affected by numerical errors 
larger than the width of the lines shown in the graph. 
The overall behavior of the approximate solution given 
by Eq. (18) is in good agreement with the numerical 
results in all cases shown. However, at low fillings, fast 
oscillations with frequency equal to ft [Fig. 8(b)] are 
superimposed onto the resonant frequency oscillations. 
These high frequencies are due to the single-particle 
BO and are not taken into account by our projection 
onto resonant families of states. However, these BO 
can be suppressed by increasing the filling [as seen in 
Fig. 8(b)] or by using larger values of f2 while ensuring 
the resonant condition. This is explained by the fact 
that the amplitude of the BO is oc At larger 

fillings, the approximation Eq. (18) breaks down since 
clusters larger than the [1111) states become relevant. 

With this, we conclude our treatment of the dynamics 
after a sudden tilt of an initial Mott-insulating state and 
turn now to the question of how to enhance transport in 
these systems. 



VI. ENGINEERING TRANSPORT: A SLINKY 
SCHEME 

In this section, we apply the projection onto resonant 
families of small cluster states to the problem of enhanc- 
ing transport of atoms on optical lattices. In particular, 
we treat systems which are prepared so that only decou- 
pled 1 11) clusters are present in the initial state. This 
can be realized using a pattern loading scheme in which 
an optical superlattice is generated by superimposing two 
lattices with different periodicity. Such a spatially selec- 
tive loading of particles onto an optical lattice has been 
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FIG. 9. (Color online) (a) Initial state at to obtained by 
a pattern loading scheme and 'slinky motion' obtained by a 
stroboscopic modulation of the lattice depth with the two 
resonant frequencies, uj = (7o ± f^o- The time sequence at 
which u is alternated is: to = 0, Jt\ = t2 = 2ti, and 
ts = 3ti. One can compare the motion of the particles to 
that of a toy slinky tumbling down a set of stairs as depicted 
above, (b) CM motion due to the stroboscopic modulation of 
the lattice depth. The different colors indicate the intervals in 
which a; = ?7o ± f^o, respectively. The dashed curve displays 
the dynamics obtained by fully diagonalizing a BHM, Eq. (1), 
with L — 7 and A'^ = 4 in the approximation U{t) = Uo, 
n{t) = Qo, and J{t) « Jsin(a;t), captured by Eq. (23). The 
solid line displays results for the same system but in which 
J{t), U{t), and fl(t) are obtained from Wannier orbitals for 
Vo,x = 0.5 and v = 0.3. The time sequence in this case is 
Jti ~ 1.43, t2 = 2ti, and ts = 3ti, where Jti is identified 
numerically as the time at which the first minimum in Xcmit) 
appears. The numbers on the top axis refer to the exact 
dynamics (solid lines). 



experimentally achieved, as discussed in Ref. 32. 

The transport through the system now is achieved by 
applying time-dependent fields. Driven tunneling by us- 
ing time-depedent fields has been addressed at the single- 



particle level before [47, 53-55]. For instance, such an 
effect has been realized in Ref. 55, where shaking the 
lattice, i.e., by applying a time-dependent linear field 
f2(i) = Oo(l -I- 7sin(wt)) and tuning the frequency close 
to that of the BO leads to an enhancement of the CM 
motion, called "super" Bloch oscillations. Similar to the 
approach of Ref. 56, we stroboscopically apply two os- 
cillatory driving fields in order to enhance the ampli- 
tude of the oscillations. At the many-body level, we 
propose to enhance transport of the atoms by perform- 
ing an amplitude modulation along the lattice direction, 
Vx{t) = Vba;(l + vsm{ujt)) where v <^ 1 and T4 the am- 
plitude of the lattice potential. The effect of the mod- 
ulation can be addressed by assuming that the Wannier 
functions possess a Gaussian profile on each lattice site. 
Using the Gaussian approximation [57], the parameters 
of the BHM are modified and one obtains [58] 

U ^ Uo [I + v sin {ujt)f\ (20) 
f2 = f2o [1 + sin (wt)] . 

The linear approximation of the Hamiltonian around 
V — results in J{t) « Jg -f Jsin(wt), where J — 

-tr T „J dlnU I dlnji ^ 

In the following, we apply the pattern loading scheme 
to realize an initial state at time to as shown in Fig. 9(a). 
In this state, we position two adjacent particles on ev- 
ery three lattice sites. For » J, BO are suppressed 
so that only resonant motion within the |11) cluster is 
possible. Hence, we can treat the dynamics to a good 
approximation by projecting onto |11) clusters. In addi- 
tion, we require Uq ^ J and |[/q — f2o| ^ J. Then the two 
energies Uo — flo and Uq + ^o are well separated from each 
other. Only by adjusting the frequency of the modulation 
to either Uq — ^Iq or Uo + flo, resonance of the |11) state 
with a particle hole excitation can be achieved; in this 
way, the |11) clusters can be resonant with either the |20) 
or the 1 02) configurations. Hence, it is possible to control 
the direction of the motion of the atoms by tuning the 
frequency accordingly, and the time evolution of the two 
particles in a cluster is either lipit)) = Cg{t)\ll)+Ce{t)\20) 
or \i){t)) = Cg{t)\n) -I- Ce(i)|02), respectively. Note that 
for a doubly occupied site, the effect is reversed, and ap- 
plying the same driving frequency will lead to a motion 
in the direction opposite to the motion induced on a pair 
of neighboring particles. In the following, we explore this 
to formulate our proposal for a transport scheme. 

In order to provide a more quantitative description of 
this behavior we first set J{t) ~ Jsm{ujt). Since we 
have only |11) clusters we can restrict the dynamics to a 
two-level system described by the effective spin model of 
Eq. (8) which in this case is 

Heff = ~V2J sin (wt) (7x + ^f^z' (21) 
where Auj — Uo±Qq. In the rotating wave approximation 
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the Hamiltonian is 



H, 



ff 



S 
2 



V2 



Jay 



(22) 



where we have set S = Acj — lo, and we assume d <C 
Aw. In this two-state representation, the center of mass 
observable is Xcm ^ Assuming an initial state in 
which Cg(0) — 1, we obtain 



'^cm (^) ^ ^ 



J2 cos (tV2P + (52 



4 J2 + 2,52 



(23) 



with the ± indicating motion up or down. Even though 
the CM amphtude is 2j2/(4j2 + 26^) < 1/2, it is never- 
theless possible to implement transport through the lat- 
tice by stroboscopically alternating the modulation fre- 
quencies Lo = Uq ± Qq in intervals of At = 7r/\^J at 
zero detuning: in this way, we transform the initial |11) 
cluster to a double occupancy, and due to the alterna- 
tion of the frequency, this is further transformed into a 
1 11) cluster shifted by one lattice spacing relative to the 
original cluster. 

Alternating w hence leads to a slinky like motion as 
depicted in Fig. 9(a), inducing a net transport of par- 
ticles through the system. Note that the net transport 
can also be uphill. In Fig. 9(b) we show as a proof of 
principle results for such a slinky motion. We display 
the exact numerical time evolution for a Bose-Hubbard 
system of L = 7 sites and = 4 particles when strobo- 
scopically modulating u and compare it to the result of an 
approximation in which U and are constant in time and 
J{t) Jsin(wt). This approximation leads to a CM mo- 
tion captured by Eq. (23). As can be seen, in both cases 
the CM motion on the time scale treated is strongly en- 
hanced. At the end of the time evolution shown, the ap- 
proximation shows transport by Aa;cm ~ 1.5 lattice spac- 
ings, while the exact solution shows Aa;cm ~ 0.9 lattice 
spacings - note that the usual CM motion is restricted to 
Axcm < 0.5 lattice spacings. Despite the difference be- 
tween the exact result and the approximate treatment. 
Fig. 9 shows that the description in terms of the slinky 
motion compares qualitatively. We therefore expect that 
for larger systems transport through the lattice should 
be realizable. In addition, we expect that for our sim- 
ple example, the transport can further be enhanced by 
optimizing the parameters. 

As mentioned in Sec. V C, we have neglected the high- 
frequency BO which will lead to a dephasing of the slinky 
motion of the atoms. This, however, can be controlled 
by choosing H. sufficiently large to dampen the BO as 
discussed in Sec. V C. Nevertheless, the tilt must be weak 
enough so that the description of the system by a one- 
band model remains valid. 



VII. SUMMARY 

We have investigated the resonant dynamics of 
strongly-interacting bosonic particles on a one- 



dimensional tilted optical lattice. At commensurate 
fillings when tuning to resonance U = ft, we find CM 
oscillations enhanced compared to the standard BO 
exhibited by non-interacting atoms. Following Ref. 20, 
the resonant dynamics can be captured by an effective 
spin- 1/2 model. Interestingly, we find signals of the crit- 
ical point in the dynamics of the CM oscillations in the 
spin model which calls for further studies. We develop 
a method to describe the dynamics at low fillings based 
on projections onto small clusters. We find that this 
approach provides a good description of the dynamics 
up to fillings n « 2/3. Using this approximation, we 
propose a scheme to engineer transport in the lattice 
by stroboscopically applying amplitude modulated 
frequencies of the lattice and envisage that this scheme 
can be realized in ongoing experiments on this system 
[21]. 

In future work, it would be interesting to explore ex- 
tensions of this work to study transport in more com- 
plex systems, such as systems with higher dimensional- 
ity [59, 60], or systems utilizing atoms that have internal 
structure [61]. In two dimensions, e.g., with a tilt along 
either one or both directions, the transport is likely to 
be quite different to that seen here in one dimension. 
Furthermore, a natural extension will be to extend these 
results on one dimensional wires to consider the trans- 
port through Atomtronic transistors [2] that naturally 
contain three ports (a base, collector, and emitter) and 
thus possess a more complicated topology. 
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Appendix A: CM motion in the fermionized regime 

In this appendix, we derive Eq. (5). We assume a lat- 
tice of infinite size centered at j = 0; the time evolution 
of an initial single particle wavefunction \ijj'^{t = 0)) then 
is 



(Al) 



n— — 00 



with the Wannier-Stark states |0„) = J2'jL-ca Jj-n{c()\j) 
and a = 2J/0. With this we obtain for the CM motion 



j 



12 



(A2) 

Using the recurrence relation, J„_i(a) + J„+i(q;) = 
^ J„(a), and the completeness relation of Bessel func- 
tions, J2j Jj-n{oi)Jj-mio:) = 5n,m, we obtain the iden- 
tity 



first order correction on C2(t) is 

ici^\t) = -V2Jc^°\t) + 2nci^\t), (B3) 
with the initial condition c^'' (0) — 0. This gives 



%/2(J2 -2^)2) 



(A3) 



With this, the CM motion of Eq. (A2) takes the form 



2n - 2e2^"^17 cos (y2Jt^ + iV2e^'*^Jsm [V^Jt^ 

(B4) 

Finally, we address the first order correction to cq {t) due 
to C2^^(t) by going back to the initial set of equations, 
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(A4) 



The CM motion for the A^-body system is obtained by 
summing over the single particle states q, resulting in 
Eq. (5). This expression is exact for an infinite system. 
However, for finite systems for ^ J wc find that this 
provides an excellent approximation even for system sizes 
as small as L = 10. Hence, the results from directly diag- 
onalizing the single particle Hamiltonian shown in Fig. 1 
and the results obtained from Eq. (5) are essentially iden- 
tical. 



Appendix B: Time-dependent perturbation theory 
for a double well system 

The time evolution of the double-well system at reso- 
nance and for strong interactions J <C is obtained by 
solving a set of coupled differential equations, 



lCo(t)=-V2[Jci{t)+C2{t)], 

ici{t) = -V2Jco{t), 

ic2{t) = -\/2Jco(i) -I- 2VLc2{t), 



(Bl) 



with initial conditions co(0) = 1 and ci(0) = C2(0) = 0. 
The non-resonant state possesses a very small population 
J^/ri^) throughout the time evolution [see Fig. 10(b)] 

which justifies the zeroth-order approximation c^\t) = 
0. The set of equations then simplifies to 



{t) = -V2Jc}-^\t) 



(B2) 



with the initial condition Cf^\o) = 1. Hence, Cg''{t) 



COS {V2Jt) and c^°^ (t) = i sin (y/2Jt) . With this, the 



(B5) 



with the initial conditions Cq^-* (0) = 0. The solution with 
the corrections gives 



Je 



co{t) = COS (V2Jt) , 
X - J + e^**" Jcos (^V2Jt^ - iV2e2'"^Osin (y2Jt^ 
Ci(t) = isin ^\/2Jt 



C2{t) 



/2(J2-2f]2) 



X 2VL- 2e^'^^^n cos (y2Jtj + iV2e^^^"jsm (^V2Jt 

(B6) 

The CM evolution then is 

3|co(0|2 + 2|ci(t)|2 + 4|c,(t)|2 

XCM{t) = 



2 (J2- 2172)2 L 

- (3 - 7J^n^ + 2n^) COS (2V2JtJ 

- 4 J2 (3 + n^) COS (V2Jt^ cos {2m) 
-UV2J^n sin (V2Jt) sin {2m) . 



(B7) 



Since we are in the regime J ^ 12, we neglect terms 
beyond to finally obtain Eq. (6), 

5 cos(2V2Jt) 



81]2 



1 + 7 



cos (^2V2 Jt^ 4cos (^V2Jtj cos {2m) 



Appendix C: Coefficients in the cluster dynamics 

We consider a homogeneous system of hard-core 
bosons at incommensurate filling {N < L). Using open 
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boundary conditions, the single-particle wavefunctions 
are ipnixi) = j^sm{^^^) where n = 1, 2, . . . , L. The 
many-body ground state then is 



iV'o) = E 



{k} 



N 

det ipn{h) 



l^) = E^wl^) (CI) 



{k} 



where the sum is over all permutations of the positions 
of the particles {k} and the determinant is 



N 

det V'n(fcj) 



V'l(fcl) -01(^2) 
-02(^1) -02(^2) 



1p2ikN) 



N 



(C2) 



For a given k, the determinant of this alternant matrix 
takes the form 



AT / 2 

det ipnih) = 



n,i—l 



L + 1 



l/N 



^_2)^(^-i)/2 



X I 1 I COS ^ — cos ^ _^ I X sm ■ 



n 



L + l 



L 



i=l 



L+l 



At low filling factors, Eq. (C3) shows that particles tend 
not to form large clusters, validating the approximation 
used in Eq. (17). The coefficients in Eqs. (17) and (18) 
are then obtained as 



FIG. 10. (Color online) Coefficients (a) |co(^)|^ \ci{t)\^ and 
(b) |c2(t)P obtained by solving Eq. (Bl) with J — 3.0 and 
Q, — 40. Note that in (b) the magnitude of |c2(t)|^ is ~ 
1% of the values shown in (a). The black solid lines depict 
the exact numerical results whereas the dashed colored lines 
depict the approximate analytic solutions. The CM motion 
in (c) compares the exact results [black solid line], results 
of Eq. (B7) which keeps terms higher than 0{J^ /Q^) [green 
dashes], and the results of Eq. (6) which keeps only terms up 
to 0{J^/n^) [red dashes]. 



AT 

det V'n(fci) 



The coefficient dI"^ defined in Sec. IV is obtained nu- 

k 

merically by taking a configuration k and counting the 
occurence of clusters \n). The n particle clusters are iden- 
tified by their contiguous empty sites. 



(C4) 
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